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The free energy and other thermodynamic properties of hexagonal-close-packed iron are calcu- 
' lated by direct ab initio methods over a wide range of pressures and temperatures relevant to the 

I Earth's core. The ab initio calculations are based on density-functional theory in the generalised- 

^\ ■ gradient approximation, and are performed using the projector augmented wave (PAW) approach. 

Thermal excitation of electrons is fully included. The Helmholtz free energy consists of three parts, 
^jQ, associated with the rigid perfect lattice, harmonic lattice vibrations, and anharmonic contributions, 

■ and the technical problems of calculating these parts to high precision are investigated. The har- 

monic part is obtained by computing the phonon frequencies over the entire Brillouin zone, and 
by summation of the free-energy contributions associated with the phonon modes. The anhar- 
' monic part is computed by the technique of thermodynamic integration using carefully designed 

' reference systems. Detailed results are presented for the pressure, specific heat, bulk modulus, ex- 

pansion coefficient and Griineisen parameter, and comparisons are made with values obtained from 
diamond-anvil-cell and shock experiments. 
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I. INTRODUCTION 



Ab-initio techniques based on density-functional theory (DPT) have played a key role for several years in the study 
of matter under extreme conditions ||l|. With recent progress in the direct ab initio calculation of thermodynamic 
free energies |^-^ , there is now great scope for the systematic and accurate calculation of thermodynamic properties 
: over a wide range of conditions. We present here extensive DFT calculations of the free energy of hexagonal-close- 
(~| packed (h.c.p.) iron under Earth's core conditions, which we have used to obtain results for a number of other 
O ■ thermodynamic quantities, including the bulk modulus, expansion coefficient, specific heat and Griineisen parameter. 
^ ' For some of these we can make direct comparisons with experimental data, which support the accuracy and realism 
^ : of the calculations; for others, the calculations supply information that is not yet available from experiments. An 
important ambition of the work is to determine thermodynamic functions without making any significant statistical- 
mechanical or electronic-structure approximations, other than those required by DFT itself, and we shall argue that 
■ we come close to achieving this. The techniques we have developed are rather general, and we believe they will find 
application to many other problems concerning matter under extreme conditions. 

The importance of understanding the high-pressure and high-temperature properties of iron can be appreciated 
by recalling that the Earth's core accounts for about 30 % of the mass of the entire Earth, and consists mainly of 
iron In fact, the liquid outer core is somewhat less dense than pure iron, and is generally accepted to contain 
light impurities such as S, O, Si or H probably the density of the solid inner core is also significantly reduced 
by impurities 0J^. Nevertheless, the thermodynamic properties of pure iron are fundamental to understanding the 
more complex real material in the core, and a large experimental effort has been devoted to measuring them. The 
difficulties are severe, because the pressure range of interest extends from 100 GPa up to nearly 400 GPa, and the 
temperature goes from ca. 3000 K to perhaps 7000 K - the temperature at the centre of the core is subject to an 
uncertainty of at least 1000 K. 

Static compression experiments with the diamond anvil cell (DAC) have been performed on Fe up to 300 GPa 
at room temperature and DAC experiments at temperatures as high as 3700 K have been reported up to 
200 GPa llo|-[l^. Our present knowledge of the phase diagram of Fe comes mainly from these experiments, though 
there are still controversies. Suffice it say that for pressures p above ca. 60 GPa and temperatures T below ca. 1500 K it 
is generally accepted that the stable phase is hexagonal close packed (h.c.p.). Recent DAC diffraction experiments jlj] 
indicate that h.c.p. is actually stable for all temperatures up to the melting line for p > 60 GPa, but earlier work 
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claimed that there is another phase, usually called (3, in a region below the melting line for pressures above ca. 40 GPa. 
The existing evidence suggests that, if the /3-phase is thermodynamically stable, its structure could either be double- 
h.c.p. (ijjl^ or orthorhombicj 17j , and in either case is closely related to the usual h.c.p. structure. According to 
very recent theoretical work [^|, h.c.p. is thermodynamically slightly more stable than double-h.c.p. at Earth's core 
pressures and temperatures. The evidence for the stability of h.c.p. over much of the high-tcmperature/high-pressure 
phase diagram is our motivation for concentrating on this phase in the present work. 

DAC measurements have given some information about thermodynamic quantities up to pressures of a few tens 
of GPa, but beyond this shock experiments (see e.g. Refs. [p0|-[2^] ) have no competitors. These experiments give 
direct values of the pressure as a function of volume |2^] on the Hugoniot curve, and have also been used to obtain 
information about the adiabatic bulk modulus and some other thermodynamic quantities on this curve. These data 
will be important in validating our calculations. Temperature is difficult to measure reliably in shock experiments [ p2| , 
and we believe that our ah initio results may be valuable in providing the needed calibration. 

The difhculties and uncertainties of experiments have stimulated many theoretical efforts. Some of the theoretical 
work has been based on simple atomistic models, such as a representation of the total energy as a sum of pair 
potentials ||2^, or the more sophisticated embedded-atom model |Q. Such models can be useful, but for accuracy 
and reliability they cannot match high-quality ah initio calculations based on density- functional theory (DFT) [p5| . 
The accuracy of DFT depends very much on the approximation used for the electronic exchange-correlation energy 
i?xc- It is known that the local density approximation (LDA) is not fully satisfactory for Fe |^6|, but that modern 
generalised-gradient approximations (GGA) reproduce a wide range of properties very accurately. These include the 
equilibrium lattice parameter, bulk modulus and magnetic moment of body-centred cubic (b.c.c.) Fe at ambient 
pressures [p7|-p9| and the phonon dispersion relations of the b.c.c. phase [|l9|. There has been much DFT work 
on different crystal structures of Fe at high pressures, and experimental low-temperature results for the pressure 
as a function of volume p{V) up to p = 300 GPa for the hexagonal close packed (h.c.p.) structure are accurately 
predicted H. Further evidence for the accuracy of DFT comes from the successful prediction of the b.c.c. to h.c.p. 



transition pressure [^|28|. With ah-initio molecular dynamics, DFT calculations can also be performed on the liquid 
state, and we have reported extensive calculations both on pure liquid Fe p9|~^ and on liquid Fe/S and Fe/0 
alloys (3|-|34). 

Recently, work has been reported |^,^ in which the thermal properties of close-packed crystalline Fe under Earth's 
core conditions are calculated using ah initio methods. In fact, the work itself was based on a tight-binding represen- 
tation of the total energy, but this was parameterised using extensive ah initio data. The authors did not attempt to 
perform the statistical-mechanical calculations exactly, but instead used the so-called 'particle in a cell' approxima- 
tion p^ , in which vibrational correlations between atoms are ignored. In spite of these limitations, the work yielded 
impressive agreement with shock data. We shall make comparisons with this work at various points in the present 
paper. 

The present DFT work is based on the GGA known as Perdew-Wang 1991 The choice of functional for 

exchange-correlation energy E^c completely determines the free energy and all other thermodynamic quantities. This 
statement is important enough to be worth elaborating. It is clear that a given i?xc exactly determines the total energy 
of the system ?7(Ri, . . . Rat) for given positions {Ri} of the atoms. But through the standard formulas of statistical 
mechanics, the function ?7(Ri, . . . Rat) exactly determines the free energy. So, provided no further electronic-stucture 
approximations are made in calculating [/(Ri, . . .Rat) from i?xc, and no statistical-mechanical approximations are 
made in obtaining F from [/(Ri, . . .Rat), then Ey^c exactly determines F. The calculation of U from E^^ has been 
discussed over many years by many authors. The work presented here is based on the projector-augmented wave 
(PAW) implementation of DFT [ ^9po|] , which is an all-electron technique similar to other standard implementations 
such as full-potential linear augmented plane waves (FLAPW) 0, as well as being very closely related to the ultrasoft 
pseudopotential (USPP) method In principle, PAW allows one to compute U with any required precision for 

a given E^c- In practical tests on Fe and other systems the technique has been shown to yield results that are 
almost indistinguishable from calculations based on FLAPW, USPP and other DFT implementations - provided the 
same i?xc is used, of course. We aim to demonstrate in this work that F can also be computed from the ah-initio 
C/(Ri, . . .Rat) to any required precision. In this sense, all the approximations made in calculating thermodynamic 
quantities are completely controlled, with the sole exception of i?xc itself. 

To clarify the precision we are aiming for in the calculation of _F, we need to explain that one of the wider 
objectives of this work is the ah initio determination of the high-pressure melting properties of Fe, to be reported 
in detail elsewhere Our approach to melting starts from the basic principle that at coexistence the Gibbs free 
energies Gso\{p,T) and Giiq(p, T) of solid and liquid are equal. But for a given pressure, the curves of Gsoiip,T) and 
G'iiq(p, r) cross at a shallow angle. The difference of slopes {dGsoi/dT)p = Sgoi and {dG\iq/dT)p — S\ic^ is equal to 
the entropy of fusion 5*111 = •S'liq — S'soi, which is comparable to /cb per atom. This means that to get the melting 
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temperature within an error of 5T, the non-canceUing errors in Gsoi and Guq must not exceed ca. k-Q^T. Ideally, we 
should like to calculate the melting temperature to within ca. 100 K, so that non-cancelling errors must be reduced 
to the level of ca. 10 meV. Our original ambition for the present work on h.c.p. Fe was to obtain F from the given 
ab-initio J7(Ri, . . . Rjv) to this precision, and to demonstrate that this has been achieved. As we shall see, this target 
has probably not been attained, but we miss it by only small factor, which will be estimated. 

We shall present results for thermodynamic quantities for pressures 50 < p < 400 GPa and temperatures 2000 < 
T < 6000 K. This is a far wider range than is strictly needed for understanding the inner core, where pressures span 
the range 330 < p < 364 GPa and T is believed to be in the region of 5000 — 6000 K. However, the wider range is 
essential in making comparisons with the available laboratory data. We set the lower limit of 2000 K for our T range 
because this is the lowest T that has been proposed for equilibrium between the h.c.p. crystal and the liquid (at lower 
T, melting occurs from the f.c.c. phase). 

In the next Section, we summarize the ab initio techniques, and give a detailed explanation of the statistical- 
mechanical techniques. The three Sections after that present our investigations of the three main components of 
the free energy, associated with the rigid perfect lattice, harmonic lattice vibrations, and anharmonic contributions, 
probing in each case the technical measures that must be taken to achieve our target precision. Sec. 6 reports our 
results for all the thermodynamic quantities derived from the free energy, with comparisons wherever possible with 
experimental measurements and previous theoretical values. Overall discussion and conclusions are given in Sec. 7. 
The implications of our results for deepening our understanding of the Earth's core will be analysed elsewhere. 



II. TECHNIQUES 
A. Ab initio techniques 

The use of DFT to calculate the energetics of many-atom systems has been extensively reviewed However, 
a special feature of the present work is that thermal electronic excitations are crucially important, and we need to 
clarify the theoretical framework in which the calculations are done. 

A fundamental quantity in this work is the electronic free energy [/(Ri, . . .Rjv;?^!) calculated at electronic tem- 
perature Tci with the N nuclei fixed at positions Ri, . . . Rat. Throughout most of the work, the statistical mechanics 
of the nuclei will be treated in the classical limit. We lose nothing by doing this, since we shall demonstrate later 
that quantum corrections are completely negligible under the conditions of interest. The Hclmholtz free energy of the 
whole system is then: 

F = -kBT\nl^j^^ J dRi...dRjv exp[-/3[/(Ri,...RAr;Tei)]| , (1) 

where A — h/{2T:MkBTy^'^ is the thermal wavelength, with M being the nuclear mass, and /3 = l/fceT. In practice, 
the electronic and nuclear degrees of freedom are in thermal equilibrium with each other, so that Td = T, but it will 
be useful to preserve the logical distinction between the two. Although t/(Ri, . . .RAr;Tci) is actually a free energy, 
we will generally refer to it simply as the total-energy function, to avoid confusion with the overall free energy F. 
It is clear from Eqn (|l]) that in the calculation of F, and hence of all other thermodynamic quantities, it makes no 
difference that [/ is a free energy; it is simply the object that plays the role of the total energy in the statistical 
mechanics of the nuclei. 

The DFT formulation of the electronic free energy U has been standard since the work of Mermin Q| , and has 
frequently been used in practical calculations [^ , ^ , though usually as a mere technical device for economising on 
Brillouin-zone sampling, rather than because electronic thermal excitation was important. The essence of finite- 
temperature DFT is that U is given by: 

U ^E~TS , (2) 

where the DFT total energy E is the usual sum of kinetic, electron- nucleus, Hartree and exchange-correlation energy 
terms, and S is the electronic entropy: 

S = -2fcB [h In/, + (1 - /,) ln(l - /O] , (3) 

i 

with fi the thermal (Fermi-Dirac) occupation number of orbital i. The electronic kinetic energy and other parts 
of E also contain the occupation numbers. We point out that in exact DFT theory the exchange-correlation (free) 
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energy E^c has an explicit temperature dependence. Very little is known about its dependence on temperature, and 
we assume throughout this work that i?xc has its zero-temperature form. 

The PAW implementation of DFT has been described in detail in previous papers . The present calculations 

were done using the VASP code | |4^ , ^. The details of the core radii, augmentation-charge cut-offs, etc are exactly as 
in our recent PAW work on liquid Fe |31|]. Our division into valence and core states is also the same: the 3p electrons 
are treated as core states, but their response to the high compression is represented by an effective pair potential, 
with the latter constructed using PAW calculations in which the 3p states are explicitly included as valence states. 
Further technical details are as follows. All the calculations are based on the form of GGA known as Perdew-Wang 
1991 p7| , p8[ . Brillouin-zone sampling was performed using Monkhorst-Pack special points and the detailed form 
of sampling will be noted where appropriate. The plane-wave cut-off of 300 eV was used, exactly as our PAW work 
on liquid Fe. 



B. Components of the free energy 

Our ab initio calculations of thermodynamic properties are based on a separation of the Helmholtz free energy F 
into the three components mentioned in the Introduction, which are associated with the rigid perfect crystal, harmonic 
lattice vibrations, and anharmonic contributions. 

To explain this separation, we start from the expression for F given in Eqn (^. We let FpGrf(7ci) = 
U(R,i, . . .R5^;Toi) denote the total free energy of the system when all atoms are fixed at their perfect-lattice po- 
sitions R5, and write C/(Ri, . . . Rat; Toi) = ^pcrf(Toi) + C^vib(Ri, • • ■ Rw; JLi), which defines the vibrational energy 
Uvih- Then it follows from eqn (|l|) that 

F = Fp,,f + i^vib , (4) 
where the vibrational free energy _Fvib is given by: 



vib 



= -A:BTln|-^ J dRi . . . dRjv exp[-/3i7vib(Ri, ■ • ■ Rjv; rci)]| • (5) 



(Note that we now omit the factor iV! from the partition function, since every atom is assumed to be confined to 
its own lattice site.) The vibrational energy Uvib can be further separated into harmonic and anharmonic parts 
(CA^ib = t^harm + C4,nharm), in tcrms of which wc Can define the harmonic vibrational free energy Fharm: 



^harm = -ksT In | J dRi . . . dRN exp[-/3C/harm(Ri , • ■ ■ Rw; Tel)] I 



(6) 



with the anharmonic free energy being the remainder Fanhaim = -Fvib — -Fharm- The harmonic energy C/harm is defined 
in the obvious way: 

?/harm = ^ I] ' (V/VjC/) • Uj , (7) 

where u/ is the displacement of atom / from its perfect-lattice position (u/ = R/ — R^) and the double gradient of 
the ab initio total energy is evaluated with all atoms at their perfect-lattice positions. Since we are dealing with a 
crystal, we shall usually prefer to rewrite J/harm in the more explicit form: 

C^harm = ^ Uisa^lsa,l'tl3Ul'tl3 , (8) 
lsa,l'tl3 

where uisa is the ath Cartesian component of the displacement of atom number s in primitive cell number Z, and 
^isa,i't/3 is the force-constant matrix. It should be noted that the present separation of F does not represent a 
separation into electronic and nuclear contributions, since thermal electronic excitations influence all three parts of 
F. 

Since all other thermodynamic functions can be obtained by taking appropriate derivatives of the Helmholtz free 
energy, the separation of F into components implies a similar separation of other quantities. For example, the pressure 
is Ppcrf + Pharm + Panharmi whcrc Ppcif — {dFpQYf/ dV)T j and similarly for the components pharm 

and Panharm- 
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C. Phonon frequencies 



The free energy of a harmonic osciUator of frequency uj is k-oTln (exp(i/3?iw) — exp(— , which has the high- 
temperature expansion k-QT h\{phLo) + k^T [^(/3?ia;)^ + 0{{fih,ujY)\ , so that the harmonic free energy per atom of 
the vibrating crystal in the classical limit is: 

Fharm = ^ ln(/3?iC^k.) , (9) 

where Wks is the frequency of phonon branch s at wavevector k and the sum goes over the first Brillouin zone, with 
N]s_s the total number of /c-points and branches in the sum. It will sometimes be useful to express this in terms of the 
geometric average Cu of the phonon frequencies, defined as: 

'^'^ = ]^El^(^ks), (10) 
ks 

which allows us to write: 

Fharm = 3fcBTln(/3to) . (11) 

The central quantity in the calculation of the frequencies is the force-constant matrix ^isa,i'tf3, since the frequencies 
at wavevector k are the eigenvalues of the dynamical matrix Dsa,tf3, defined as: 

DscM^) = ]g E • (^"* - ^is)] , (12) 

where Rj*^ is the perfect-lattice position of atom s in primitive cell number I. If we have the complete force-constant 
matrix, then Dsa,ti3 and hence the frequencies Wks can be obtained at any k, so that Co can be computed to any 
required precision. In principle, the elements of ^isa,i'ti3 are non-zero for arbitrarily large separations | R^^ — R^^ |, 
but in practice they decay rapidly with separation, so that a key issue in achieving our target precision is the cut-off 
distance beyond which the elements can be neglected. 

We calculate ^isa,i'tf3 by the small-displacement method, in a way similar to that described in Ref. ||5^. The basic 
principle is that ^isa,i'ti3 describes the proportionality between displacements and forces. If the atoms are given small 
displacements uisa from their perfect-lattice positions, then to linear order the forces Fisa are: 

Flsa ^ -^^^lsa,l'tl3Ul'tf3 ■ (13) 

Within the ab initio scheme, all the elements ^isa.i'tp are obtained for a given I'tp by introducing a small displacement 
ui'tp, all other displacements being zero, minimising the electronic free energy, and evaluating all the forces Fisa- In 
practice, the displacement amplitude ui'tp must be made small enough to ensure linearity to the required precision, 
and this sets the precision with which the electronic free energy must be minimised. 

By translational symmetry, the entire force-constant matrix is obtained by making three independent displacements 
for each atom in the primitive cell, and this means that no more than SiVbas calculations are needed, where iVbas is 
the number of atoms in the primitive cell. This number can be reduced by symmetry. If, as in the h.c.p. crystal, 
all atoms in the primitive cell are equivalent under operations of the space group, then the entire force-constant 
matrix can be obtained by making at most three displacements of a single atom in the primitive cell: from ^isa,i'tp 
for one chosen atom I't, one obtains ^isa,i'tf3 for all other I't. Point-group symmetry reduces the number still 
further if linearly independent displacements of the chosen atom are equivalent by symmetry. This is the case in 
the h.c.p. structure, since displacements in the basal plane related by rotations about the c-axis by ±120° are 
equivalent by symmetry; this means that two calculations, one with the displacement along the c-axis, and other 
with the displacement in the basal plane, suffice to obtain the entire ^isa.i'tfs matrix. The basal-plane displacement 
should be made along a symmetry direction, because the symmetry makes the calculations more efficient. Since the 
exact $;sQ,;'t/3 matrix has point-group symmetries, the calculated ^isa.i'tp must be symmetrised to ensure that these 
symmetries are respected. The symmetrisation also serves to eliminate the lowest-order non-linearities in the relation 
between forces and displacements Q . 

It is important to appreciate that the ^isa.i'tf} in the formula for Dsa.tpi^) is the force-constant matrix in the 
infinite lattice, with no restriction on the wavevector k, whereas the ab initio calculations of $;sa.i't/3 can only be done 
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in supercell geometry. Without a further assumption, it is strictly impossible to extract the infinite-lattice ^isa,i'tf3 
from supercell calculations, since the latter deliver information only at wavevectors that are reciprocal lattice vectors 
of the superlattice. The further assumption needed is that the infinite-lattice ^isa,i'ti3 vanishes when the separation 
Ri't — Ris is such that the positions Rjs and Ri't lie in different Wigner-Seitz (WS) cells of the chosen superlattice. 
More precisely, if we take the WS cell centred on Hin, then the infinite-lattice value of ^isa,i't0 vanishes if R;^ is 
in a different WS cell; it is equal to the supercell value if R/s is wholly within the same WS cell; and it is equal to 
the supercell value divided by an integer P if R;^ lies on the boundary of the same WS cell, where P is the number 
of WS cells having R/s on their boundary. With this assumption, the ^isa,vtp elements will converge to the correct 
infinite-lattice values as the dimensions of the supercell are systematically increased. 



D. Anharmonicity 

1. Thermodynamic integration 

Although we shall show that the anharmonic free energy i^anhaim is numerically fairly small, it is far more challenging 
to calculate than i^perf or i^harm, because there is no simple formula like eqn (|9|), and the direct computation of the 
multi-dimensional integrals in the free-energy formulas such as eqn is impossible. Instead, we use the technique of 
thermodynamic integration (see e.g. Ref. to obtain the difference i^vib— ^harm, as developed in earlier papers 

Thermodynamic integration is a completely general technique for determining the difference of free energies Fi — Fq 
for two systems whose total-energy functions are Ui and Uq. The basic idea is that Fi — Fq represents the reversible 
work done on continuously and isothermally switching the energy function from Uq to Ui. To do this switching, a 
continuously variable energy function U\ is defined as: 

Ux = il- X)Uo + XUi , (14) 

so that the energy goes from C/q to C/i as A goes from to 1. In classical statistical mechanics, the work done in an 
infinitesimal change d\ is: 

dF = {dUx/dX)xdX ^ {Ui - Uo)xdX , (15) 
where {■)x represents the thermal average evaluated for the system governed by Ux- It follows that: 

Fi-Fo= [ dX{Ui-Uo)x. (16) 
Jo 

In practice, this formula can be applied by calculating {Ui — Uo)x for a suitable set of A values and performing the 
integration numerically. The average (Ui — Uo)x is evaluated by sampling over configuration space. 

For the anharmonic free energy, a possible approach is to choose Uq as C/harm and Ui as Uvih , so that Fi — Fq is the 
anharmonic free energy i^anharm- This was the procedure used in our earlier ab initio work on the melting of Al 
and a related technique was used by Sugino and Car in their work on Si melting. However, the calculations are 
rather heavy, and the need for extensive sampling over the electronic Brillouin zone in the ab initio calculations makes 
it difficult to achieve high precision. We have now developed a more efficient two-step procedure, in which we go first 
from the harmonic ab initio system J/harm to an intermediate reference system U^cf which closely mimics the full ab 
initio total energy Uvib', in the second step, we go from C/ief to C/vib- The anharmonic free energy is thus represented 
as: 

-^anharm — (-^vib -^rcf ) (-^rcf -^harm) ; (-^'^) 

and the two differences are calculated by separate thermodynamic integrations: 

-Fvib — ^rcf — dX{Uvib — C^rcf)r 

Jo 

Frcf — Fharm — dX {U^-cf — UharmYx ■ (18) 

Jo 

To distinguish clearly between these two parts of the calculation, we denote by ( • )? ^^"^ thermal average taken in the 
ensemble generated by the switched total energy U^^ = (1 ^ A)C/haim + XUid and by ( ■ the corresponding average 
for C/;^ = (1 - A)C/ref + AC/vib. 
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The crucial point of this is that f/ref is required to consist of an empirical model potential which quite accurately 
represents both the harmonic and anharmonic parts of the ab initio total energy Uvih- Since it is a model potential, 
the thermodynamic integration for F^et — -Fkarm can be performed with high precision on large systems. The difference 
Fvih ~ -frefj by Contrast, involves heavy ab initio calculations, but provided a good U^cf can be found these are 
manageable. We return to the problem of searching for a good Uid in Sec. IID 3 



2. Calculation of thermal averages 

The calculation of thermal averages is just the standard problem of computational statistical mechanics, and can be 
accomplished by any method that allows us to draw unbiased samples of configurations from the appropriate ensemble. 
In this work, we employ molecular dynamics simulation. This means, for example, that to calculate (?7rcf — C^harm)A^ 
we generate a trajectory of the system using equations of motion derived from the total energy function U^^. In the 
usual way, an initial part of the trajectory is discarded for equilibration, and the remainder is used to estimate the 
average. The duration of this remainder must suffice to deliver enough independent samples to achieve the required 
statistical precision. 

The key technical problem in calculating thermal averages in nearly harmonic systems is that of ergodicity. In the 
dynamical evolution of a perfectly harmonic system, energy is never shared between different vibrational modes, so 
that a system starting at any point in phase space fails to explore the whole of phase space. This means that in a nearly 
harmonic system exploration will be very slow and inefficient, and it is difficult to generate statistically independent 
samples. We solve this following Ref. ||]: the statistical sampling is performed using Andersen molecular dynamics [ p2| , 
in which the atomic velocities are periodically randomised by drawing them from a Maxwellian distribution. This 
type of simulation generates the canonical ensemble and completely overcomes the ergodicity problem. 



3. Reference system 



The computational effort needed to calculate Fvib — -Frof is greatly reduced if the difference of total energies C^vib — f^rof 
is small, for two reasons. First, the amount of sampling needed to calculate (C/vib — Urcf)x to a given precision is 
reduced if the fluctuations of C/vib — C/rcf are small; second, the variation of (C/vib — Uicf)\ as A goes from to 1 is 
reduced. In fact, if the fluctuations are small enough, this variation can be neglected, and it is accurate enough to 
approximate _Fvib — ^rcf — (C/vib — f^ref)rof, with the average taken in the reference ensemble. If this is not good 
enough, the next approximation is readily shown to be: 

Fvib - F,,f ~ (C/vib - C/rcf)rcf - T^K^ ([U.ih " C/^cf " (C/vib " C/,cf )rcf )] ' ) ■ (19) 

ZK^l \ I rcf 

Our task is therefore to search for a model C/icf for which the fluctuations of C/vib ~ U^d are as small as possible. 

The question of reference systems for Fe has already been discussed in our recent ab initio simulation work on the 
high-pressure liquid . We showed there that a remarkably good reference model is provided by a system interacting 
through inverse-power pair potentials: 

t^iP = ^ E -^(1 - I) ' (20) 

where 0(r) — B/r°'^ with B and a adjusted to minimise the fluctuations of the difference between C/jp and the ab 
initio energy. Unfortunately, we shall show that this is an unsatisfactory reference model for the solid, because the 
harmonic phonon dispersion relations produced by C/ip differ markedly from the ab initio ones. It is a particularly 
poor reference model at low temperatures where anharmonic corrections are small, because in that regime a good 
reference system must closely resemble C/harm- However, we find that C/jp becomes an increasingly good reference 
system as T approaches the melting temperature. We therefore adopt as a general form for the reference system a 
linear combination of C/harm and C/ip : 

C/ref = CiC/harm + C2C/1P . (21) 

The coefficients ci and C2 are adjusted to minimise the intensity of the fluctuations of C/vib — C/icf for each thermody- 
namic state. 
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Now consider in more detail how this optimisation of Uret is to be done. In principle, the ensemble in which 
we have to sample the fluctuations of ?7vib — C^ref is the one generated by the continuously switched total energy 
(1 — A)C/ref + At/vib that governs the thermodynamic integration from t/rcf to [/vib- In practice, this is essentially the 
same as sampling in either of the ensembles associated with U^ei or C/vib, provided the fluctuations of t/vib — C^ref are 
indeed small. But even this poses a problem. We are reluctant to sample in the ensemble of J7vib, because extensive 
(and expensive) ab initio calculations are needed to achieve adequate statistical accuracy. On the other hand, we 
cannot sample in the ensemble of J/ref without knowing ?7ref, which is what we are trying to find. We resolve this 
problem by constructing an initial optimised t/ref by minimising the fluctuations in the ensemble of J7harm- We then 
use this initial J/rcf to generate a new set of samples, which is then used to reoptimise C/ref- In principle, we should 
probably repeat this procedure until C/ref ceases to vary, but in practice we stop after the second iteration. Note 
that even this approach requires fully converged ab initio calculations for a large set of configurations. But since the 
configurations are generated with the potential model t/ref , statistically independent samples are generated with much 
less effort than if we were using C/vib to generate them. 

III. THE RIGID PERFECT LATTICE 

The energy and pressure as functions of volume for h.c.p. Fe at low temperatures (i.e temperatures at which lattice 
vibrations and electronic excitations can be neglected) have been extensively investigated both by DAC experiments |^ 
and by DFT studies |||,|8|, including our own earlier USPP H and PAW ^ calculations. The various DPT cal- 
culations agree very closely with each other, and reproduce the experimental p{V) relation very accurately, especially 
at high pressures: the difference between our PAW pressures and the experimental values ranges from 4.5 % at 
100 GPa to 2.5 % at 300 GPa, these deviations being only slightly greater than the scatter on the experimental values. 

The present DFT calculations on the rigid perfect lattice give Fpcrf for any chosen volume and electronic temperature 
Tel. In order to achieve our target precision of 10 meV/atom for the free energy, careful attention must be paid to 
electronic Brillouin-zone sampling. All the calculations presented in this Section employ the 15 x 15 x 9 Monkhorst- 
Pack set, which gives 135 fc-points in the irreducible wedge of the Brillouin zone. The fc-point sampling errors with 
this set have been assessed by repeating selected calculations with 520 fc-points in the irreducible wedge. Tests at 
the atomic volume V = 8.67 A'^ show that the errors are 4.0, 2.0 and 0.5 meV/atom at temperatures of 500, 1000 
and 2000 respectively, so that, as expected, the errors decrease rapidly with increasing T^. Since temperatures below 
2000 K are not of interest here, we conclude that with our chosen Monkhorst-Pack set fc-point errors are completely 
negligible. 

We have done direct DFT calculations at a closely spaced set of atomic volumes going from 6.2 to 11.4 A"^ at intervals 
of 0.2 A'^, and at each of these volumes we have made calculations at Toi values going from 200 to 10,000 K at intervals of 
200 K. At every one of these state points, we obtain the value of i^pcrf • The calculations also deliver directly the internal 
energy -Epcrf and the entropy S'porf (see Eqn ^ . The specific heat is then obtained either by numerical differentiation 
of the iJpcrf results (Cpcrf — {dEpcrf/dT)v) or analytically from the formula for S'porf (Cpcrf = T{dSpcYf/dT)Y). If 
we ignore the temperature dependence of the Kohn-Sham energies, dSpcri/dT can be evaluated analytically using the 
text-book formula for the Fermi-Dirac occupations numbers: fi = l/[exp(/3(ei — m)) + 1], where is the Kohn-sham 
energy of orbital i and fj, is the electronic chemical potential. It was pointed out by Wasserman et al. p5| that the 
neglect of the dependence of on Ttoi is a very accurate approximation, and out tests confirm that the errors incurred 
by doing this are negligible. 

Our analytically calculated results for Cporf are reported in Fig. |l| for the temperature range — 6000 K at the 
atomic volumes V — 7.0, 8.0, 9.0 and 10.0 A^/atom. The key point to note is that Cporf becomes large at high 
temperatures, its value of ca. 2fcB at 6000 K being comparable with the Dulong-Petit specific heat of lattice vibrations 
(Sfce). This point about the large magnitude of the electronic specific heat has been emphasised in several previous 
papers |^,^,^ , and the values reported here are close to those calculated by Wasserman et al. , though the latter 
actually refer to f.c.c. Fe. This means that full inclusion of thermal electronic excitations, as done here, is crucial to 
a correct description of the thermodynamics of Fe at core conditions. 

The linear dependence of Cporf on T evident in Fig. |^ at low T (Cporf — "fT + 0{T^)) is what we expect from the 
standard Sommerfeld expansion |Q for electronic specific heat in powers of T, which shows that the low-temperature 
slope is given by 7 = ^■7T^k^g{Ep), where giEp) is the electronic density of states (DOS, i.e. the number of states 
per unit energy per atom) at the Fermi energy Ep. Our calculated DOS at the atomic volumes V = 7.0 and 10.0 A^ 
(Fig. 1^) shows, as expected, that the width of the electronic d-band increases on compression, so that g{Ep) and 
hence 7 decrease with decreasing atomic volume. As a cross-check, we have calculated 7 directly from the density of 
states, and we recover almost exactly the low temperature slope of Cporf- 
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In order to obtain other thermodynamic functions, we need a fit to our Fperf results. At each temperature, we fit 
the results to the standard Birch-Murnaghan form, using exactly the procedure followed in our recent work on the 
Fe/0 system Q. This involves fitting the 22 values of Fpcrf at a given temperature using four fitting parameters 
{Eq, Vq, K and K' in the notation of Ref. [0). We find that at all temperatures the r.m.s. fitting errors are less 
than 1 meV at all points. The temperature variation of the fitting parameters is then represented using a polynomial 
of sixth degree. 

Electronic excitations have a significant effect on the pressure, as can be seen by examining the T-dependence of the 
perfect-lattice pressure ppcrf — —{dFpcrildV)T- We display in Fig. || the thermal part Apporf ofpporf, i-e. the difference 
between Pperf at a given T and its zero-temperature value. The thermal excitation of electrons produces a positive 
pressure. This is what intuition would suggest, but it is worth noting the reason. Since Fpcrf (Td) = Fpcrf (0) — ^T^7(y) 
at low temperatures, the change of pressure due to electronic excitations is Ap = ^T^dj/dV in this region. But 
dj/dV > 0, so that the electronic thermal pressure must be positive. To put the magnitude of this pressure in 
context, we recall that at the Earth's inner-core boundary (ICB) the pressure is 330 GPa and the temperature is 
believed to be in the range 5000 — 6000 K, the atomic volume of Fe under these conditions being ca. 7 A'^. Our results 
then imply that electronic thermal pressure accounts for ca. 4 % of the total pressure, which is small but significant. 

IV. THE HARMONIC CRYSTAL 
A. Convergence tests 

We have undertaken extensive tests to show that our target precision of 10 meV/atom can be attained for the 
harmonic free energy Fharm- It is useful to note that at T = 6000 K an error of 10 meV represents 2 % of fceT, and 
since Fharm is given in terms of the geometric mean frequency u) by 3kBT hi{hu} / kBT) , we must achieve a precision of 
0.7 % in cj. A sufficient condition for this is that we obtain the phonon frequencies Ws(k) for all wavevectors k and 
branches s to this precision, but this may not be necessary, since there can be cancellation of errors at different k 
and/or s. Convergence of Co must be ensured with respect to four main parameters: the atom displacement used to 
calculate the force-constant matrix ^isa,i'tf3] the electronic fc-point sampling; the size of repeating cell used to obtain 
^isa,i'ti3] and the density of /c-point mesh used in calculating 6) from the by integration over the phonon Brillouin 
zone (see Eqn (||)). Convergence can be tested separately with respect to these four parameters. 

Integration over the phonon Brillouin zone is performed using Monkhorst-Pack fc-points |^^. We find that the 
set having 364 fc-points in the irreducible wedge achieves a precision in -Fharm of better than 1 meV/atom at all 
the temperatures of interest, and this fc-point set was used in the calculations presented here. The effect of atomic 
displacement amplitude was tested with the force-constant matrix generated using a 2 x 2 x 2 repeating cell at the 
atomic volume 8.67 A"^, and amplitudes ranging from 0.0148 to 0.118 A were used. The systematic variation of the 
resulting Fharm showed that with an amplitude of 0.0148 A the error is less than 1 meV/atom. 

The errors associated with electronic fc-point sampling in the calculation of the force-constant matrix were initially 
assessed with ^isa,i'ti3 obtained from the 2x2x2 repeating cell at the atomic volume V = 8.67 A'^. We found that at 
Toi = 4300 K convergence of Fharm is obtained within 2 meV/atom if the 3x3x2 Monkhorst-Pack set of electronic 
fc-points is used. This is close to being satisfactory, but starts to produce significant errors at lower temperatures. 
Our definitive calculations were actually done with the more extensive 5x5x5 Monkhorst-Pack electronic set. With 
this set, our tests performed with $/sa,i't/3 obtained from the 3x3x2 repeating cell and V — 9.17 A'^ show that the 
electronic fc-point error is now ca. 0.1 meV/atom even at Td = 1500 K. At higher electronic temperatures and with 
larger repeating cells, the error will of course be even smaller. 

Finally, we have tested the convergence of -Fharm with respect to the size of repeating cell used to generate ^isa,i'tf3- 
A wide range of different cell sizes and shapes were studied, including 2x2x2, 3x3x2,4x4x4 and 5x5x3, the 
largest of these containing 150 atoms in the repeating cell. The tests showed that with the repeating cell 3x3x2 
the error in Fharm calculated at = 8.67 A'^ at 4300 K is a little over 2 meV/atom, and we adopted this cell size for 
all our calculations of ^isa,i'ti3- We expect the error to be similar at other volumes, and to be roughly proportional 
to temperature, so that it should be insignificant over the whole range of states of interest. 

B. Dispersion relations, average frequency, free energy 

In Fig. ^ we present the harmonic phonon dispersion relations at the two atomic volumes 8.67 and 6.97 A^^ calculated 
with Tel = 4000 K. We are not aware of previous direct ab initio calculations of the phonon frequencies of high-pressure 
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h.c.p. Fe, but there are published dispersion relations derived from a 'generalised pseudopotential' parameterisation 
of FP-LMTO calculations performed by Soderlind et al. using the LDA at the atomic volume 6.82 A^. The 
agreement of their phonon frequencies with ours is far from perfect. For example, we find that the maximum frequency 
in the Brillouin zone calculated at V = 6.82 A'^ is at the F-point and is 21.2 THz, whereas they find the maximum 
frequency at the M-point with the value 17.2 THz. This is not unexpected, since they report that the generalised 
pseudopotential scheme fails to reproduce accurately some phonon frequencies calculated directly with FP-LMTO in 
the f.c.c. Fe crystal |2^; in addition, the LDA used by them is known to underestimate phonon frequencies in Fe [ p9[ . 

Casual inspection suggests that our dispersion curves at the two atomic volumes are almost identical apart from an 
overall scale factor. This suggestion can be judged from the right-hand panel of Fig. ^ where we plot as dashed curves 
the dispersion curves at V = 8.67 A'^ scaled by the factor 1.409 - the reason for choosing this factor will be explained 
below. The comparison shows that the curves at the two volumes are indeed related by a single scaling factor to 
within ca. 5 %. We also take the opportunity here to check how well the inverse-power potential model Uip (see 
Eqn (po|)) reproduces phonon frequencies. To do this, we take exactly the same parameters B and a specifying (f){r) 
that reproduced well the properties of the liquid [0, namely a = 5.86 and B such that for r = 2.0 A (f>{r) — 1.95 eV. 
The phonons calculated from this model are compared with the ab initio phonons at atomic volume V = 8.67 A in 
the left panel of Fig. ^. Although the general form of the dispersion curves is correctly reproduced, it is clear that 
the model gives only a very rough description, with discrepancies of as much as 30 % for some frequencies. 

We performed direct ab initio calculations of the dispersion relations and hence the geometric mean frequency ui 
for seven volumes spaced roughly equally from 9.72 to 6.39 A^, and for each of these volumes for T^i from 1000 
to 10,000 K at intervals of 500 K. The results for Cu as function of volume are reported in Fig. ^ for the three 
temperatures T^i — 2000, 4000 and 6000 K. We use a (natural) log-log plot to display the results, so that the 
negative slope 7ph = —dhiQ/dhiV is the so-called phonon Griineisen parameter. (The relation between 7th and the 
thermodynamic Griineisen parameter 7 will be discussed in Sec. VI B.) We note that if phonon dispersion curves at 
two different volumes are related by a simple scaling factor, this must be the ratio of Q values at the two volumes. 
The scaling factor used in Fig. ^ was obtained in exactly this way from our uj results. The Griineisen parameter 7ph 
increases with increasing volume, in accord with a widely used rule of thumb [^^. We find that 7ph goes from 1.34 
at y = 6.7 A^ to 1.70 at y = 8.3 A^, but then decreases slightly to 1.62 at F = 9.5 A^. Fig. | also allows us to 
judge the effect of Td on phonon frequencies: for all volumes studied, the frequencies decrease by ca. 4 % as T^i goes 
from 2000 to 6000 K. However, we mention that for the higher volumes, though not for the smaller ones, Co slightly 
increases again as Td goes to still higher values. To enable the ui data to be used in thermodynamic calculations, we 
parameterise the temperature dependence of IntD at each volume as a -f bT^ + cT^ + eT^, and the volume dependence 
of the four coefhcients a, 6, c and e as a third-degree polynomial in V . 

We now return to the matter of quantum nuclear corrections. Since the leading high-temperature correction to 
the free energy is ^kBT{(3huj)'^ per mode and there are three modes per atom, the quantum correction to -Fharm is 
^kBT{ph{uj'^y^'^)'^ per atom, where (w^) denotes the average of over wavevectors and branches. At the lowest 
volume of interest, V = 7 A^, (ti>^)^/^/27r is roughly 15 THz. At the lowest temperature of interest, T = 2000 K, this 
gives a quantum correction of 3 meV/atom, which is small compared with our target precision. 



C. Harmonic phonon specific heat and thermal pressure 

If the mean frequency a) were independent of temperature, the constant-volume specific heat Charm due to harmonic 
phonons would be exactly 3fcB per atom in the classical limit employed here. We find that its temperature dependence 
yields a slight increase of Charm above this value, but this is never greater than 0.25fcB under the conditions of interest. 
The harmonic phonon pressure Pharm as a function of atomic volume at different temperatures is reported in Fig. ||. 
Comparison with Fig. |^ shows that Pharm is always much bigger (by a factor of at least three) than the electronic 
thermal pressure under the conditions of interest. At ICB conditions {p = 330 GPa, T 5000 — 6000 K), pharm 
account for ca. 15 % of the total pressure. 



V. ANHARMONIC FREE ENERGY 



A. Optimisation of reference system 



It was stressed in Sec. II D 3 that optimisation of the reference system greatly improves the efficiency of the anhar- 
monic calculations. We investigated the construction of the reference system in detail at the atomic volume 8.67 A'^, 
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with the optimisations performed for a simulated system of 16 atoms. The calculation of the anharmonic free energy 
itself for a system as small as this would not be adequate, but we expect this system size to suffice for the optimisation 
of f/ref- The initial sample of configurations (see Sec. II D 3) was taken from a simulation of duration 100 ps performed 
with the total energy [/harm, with velocity randomisation typically every 0.2 ps. Configurations were taken every 1 ps, 
so that we obtain a sample of 100 configurations. In computing the energy difference Uvih — U^ei for these configu- 
rations, the ab initio energy C/vib was always computed using 5x5x3 Monkhorst-Pack electronic fc-point sampling 
(38 /c-points in the full Brillouin zone). Once the preliminary optimisation had been performed with configurations 
generated like this, the resulting C/rof was used to produce a new set of 100 configurations with an Andersen MD 
simulation of the same duration as before, and the reference system was reoptimised. 

This entire procedure was carried out at temperatures of 1000 and 4000 K. The values of the optimisation coefficients 
(see Eqn (|2l|)) were ci = 0.2, C2 = 0.8 at the high temperature and ci = 0.7, C2 = 0.3 at the low temperature. (We 
do not require that ci + C2 = 1, though this happens to be the case here.) As expected, Urd resembles C/harm quite 
closely at the low temperature and Uip quite closely at the high temperature. 

In view of the labour involved in the optimisation, we wanted to find out whether the detailed choice of ci and C2 
makes a large difference to the strength of the fluctuations of Uvih ~ C/rof- To do this, we computed these fluctuations 
at several temperatures, using the two reference models just described, i.e. without optimising the Ci coefficients at 
each temperature. Our conclusion is that the values ci = 0.2, C2 = 0.8 can safely be used at all the state points of 
interest, without incurring large fluctuations, and we therefore used this way of making the reference system in all 
subsequent calculations. 



B. From harmonic ab initio to reference to full ab initio 



The thermodynamic integration from ab initio harmonic to reference was done with nine equally-spaced A-points 
using Simpson's rule, which gives an integration precision well in excess of our target. To investigate the influence 
of system size, integration from C/harm to C/rof was performed for systems of 12 different sizes, going from 16 to 1200 
atoms. The calculations were also repeated with the force-constant matrix in C/harm generated with cells containing 
from 16 to 150 atoms (see Sec. II C). These tests showed that if the thermodynamic integration is done with a system 
of 288 atoms and the force constant used for C/harm is generated with the 36-atom cell, then the resulting difference 
-Fief ~ Fiiarm IS Converged to better than 3 meV/atom. 

To compute the difference F^ih — -Frof, we used the second-order expansion formula given in Eqn (|l|). Given the 
small size of the fluctuations of C/vib — C/ref , we expect this to be very accurate. The calculations of FVib — Flcf were all 
done with the 16-atom system. Tests with 36- and 64-atom systems show that this free-energy difference is converged 
with respect to size effects to within ca. 2 meV. 

A summary of all our results for Fj-cf — Fharm and FVib — Frcf , and the resulting values of Fanharm = Fvib — Fharm 
are reported in Table |. The anharmonic free energy is always negative, so that anharmonicity stabilises the solid. As 
expected, the anharmonic free energy is small (less than or comparable with our target precision of 10 meV/atom) 
at low temperatures, but increases rapidly at high temperatures. The temperature at which it becomes appreciable 
is higher for smaller atomic volumes. 

In classical statistical mechanics, -Fanharm is expected to go as at low temperatures, and in fact we flnd that 
Fknharm = a(V)T^ givcs a good representation of our results for all the temperatures studied. The volume dependence 
of a(y) is adequately represented by a{V) = ai + a-iV , with ol\ = 2.2 x 10^^ eV K^^ and a-i — -6.0 x 10^^" eV 
per atom. 



C. Anharmonic specific heat and pressure 

Within the parameterisation just described, the anharmonic contribution to the constant-volume speciflc heat 
Canharm is proportioual to T and varies linearly with V. As an indication of its general size, we note that Canharm 
increases from 0.09 to 0.18 /cb at 2000 K and from 0.28 to 0.53 /cb at 6000 K as F goes from 7 to 10 A^. The 
anharmonic contribution to the pressure is independent of volume, and is proportional to T^. It increases from 0.4 
to 3.5 GPa as T goes from 2000 to 6000 K, so that even at high temperatures it is barely significant. 
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VI. THERMODYNAMICS OF THE SOLID 



We now combine the parameterised forms for -Fporf, -Fharm and -Fanharm presented in the previous three Sections 
to obtain the total free energy of the h.c.p. crystal, and hence, by taking appropriate derivatives, a range of other 
thermodynamic functions, starting with those measured in shock experiments. 



A. Thermodynamics on the Hugoniot 

In a shock experiment, conservation of mass, momentum and energy require that the pressure pn, the molar internal 



energy En and the molar volume Vh in the compression wave are related by the Rankine-Hugoniot formula |57|: 

^pM - 1/h) =En-Eo, (22) 

where Eq and Vq are the internal energy and volume in the zero-pressure state before the arrival of the wave. The 
quantities directly measured are the shock-wave and material velocities, which allow the values of pn and Vh to 
be deduced. From a series of experiments, pn as a fmiction of Vh (the so-called Hugoniot) can be derived. The 
measurement of temperature in shock experiments is attempted but problematic [ p2| . 

The Hugoniot curve PhCI^h) is straightforward to compute from our results: for a given Vh, one seeks the temperature 
at which the Rankine-Hugoniot relation is satisfied; from this, one obtains pn (and, if required. En). In experiments 
on Fe, Vq and Eq refer to the zero-pressure b.c.c. crystal, and we obtain their values directly from GGA calculations, 
using exactly the same PAW technique and GGA as in the rest of the calculations. Since b.c.c. Fe is ferromagnetic, 
spin polarisation must be included, and this is treated by spin interpolation of the correlation energy due to Vosko et 
al., as described in Refs. [ pT|j40) . The value of Eq includes the harmonic vibrational energy at 300 K, calculated from 
ah initio phonon dispersion relations for ferromagnetic b.c.c. Fe. 

Our ab initio Hugoniot is compared with the measurements of Brown and McQueen pl| in Fig. |^. The agreement is 
good, with discrepancies ranging from 10 GPa at V — 7.8 A'^ to 12 GPa at V^ = 8.6 A'^ These discrepancies are only 



slightly greater than those found for the room-temperature static p{V) curve (see Sec. HI), which can be regarded 
as giving an indication of the intrinsic accuracy of the GGA itself. Another way of looking at the accuracy to be 
expected of the GGA is to recalculate the Hugoniot using the experimental value of the b.c.c. Vq (11.8 A'^, compared 
with the ab initio value of 11.55 A'^). The Hugoniot calculated in this way is also plotted in Fig. 0, and we see that 
this gives almost perfect agreement with the experimental data in the pressure range 100 — 240 GPa. We deduce from 
this that the ab initio Hugoniot deviates from the experimental data by an amount which should be expected from 
the known inaccuracies of the GGA applied to Fe. A similar comparison with the experimental Hugoniot was given 
in the tight-binding total-energy work of Wasserman et al. |35|] , and their agreement was as good as ours. We discuss 
the significance of this later. 

Our Hugoniot temperature as a function of pressure is compared with the experimental results of Brown and 
McQueen and of Yoo et al. in Fig. ||. The ab initio temperatures agree well with those of Brown and 
McQueen, but fall substantially below those of Yoo et al., and this supports the suggestion of Ref. |^ that the Yoo 
et al. measurements overestimate the Hugoniot temperature by ca. 1000 K. 

A further quantity that can be extracted from shock experiments is the bulk sound velocity as a function of atomic 
volume on the Hugoniot, which is given by ub = {Kg/ pY^^ , with Ks = —V{dp/dV)s the adiabatic bulk modulus 
and p the mass density. Since Ks can be calculated from our ab initio pressure and entropy as functions of V and T, 
our calculated Ks can be directly compared with experimental values (Fig. ^). Here, there is a greater discrepancy 
than one would wish, with the theoretical values falling significantly above the Ks values of both Refs ||2^ and [ pO[ , 
although we note that the two sets of experimental results disagree by an amount comparable with the discrepancy 
between theory and experiment. 

For what it is worth, we show in Fig. ^ a comparison between our calculated thermal expansivity on the Hugoniot 
with values extracted from shock data by Jeanloz |2^. The latter are very scattered, but is clear that the theoretical 
values have similar magnitude. However, our values vary little along the Hugoniot, whereas the experimental values 
seem to decrease rather rapidly with increasing pressure. 
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B. Other thermodynamic quantities 



We conclude our presentation of results by reporting our ab initio predictions of quantities which conveniently 
characterise h.c.p. Fe at high pressures and temperatures, and allow some further comparisons with the predictions 
of Refs 1^ and Our results are presented as a function of pressure on isotherms at T = 2000, 4000 and 6000 K. 
At each temperature, we give results only for the pressure range where, according to our ab initio melting curve, the 
h.c.p. phase is thermodynamically stable. 

The total constant-volume specific heat per atom C-u (Fig. ^l]) emphasises again the importance of electronic 
excitations. In a purely harmonic system, Cy would be equal to 3/cb, and it is striking that Cy is considerably greater 
than that even at the modest temperature of 2000 K, while at 6000 K it is nearly doubled. The decrease of Cy with 
increasing pressure evident in Fig. |ll| comes from the suppression of electronic excitations by high compression, and 
to a smaller extent from the suppression of anharmonicity. 

The thermal expansivity a (Fig. |l^) is one of the few cases where we can compare with DAC measurements [ p^ . 
The latter show that a decreases strongly with increasing pressure and the ab initio results fully confirm this. Our 
results also show that a increases significantly with temperature. Both trends are also shown by the tight-binding 
calculations of Ref. IBS], though the latter differ from ours in showing considerably larger values of a at low pressures. 
We note that Ref. |135| reported results for a at temperatures only up to 2000 K, so a full comparison is not possible. 

The product aKT of expansivity and isothermal bulk modulus is important because it is sometimes assumed to be 
independent of pressure and temperature over a wide range of conditions, and this constancy is used to extrapolate 
experimental data. Our predicted isotherms for aKx (Fig. |l^) indicate that its dependence on p is indeed weak, 
especially at low temperatures, but that its dependence on T certainly cannot be ignored, since it increases by at 
least 30 % as T goes from 2000 to 6000 K at high pressures. Wasserman et al. come to qualitatively similar 
conclusions, and they also find values of ca. 10 MPa K"'^ at T ~ 2000 K. However, it is disturbing to note that the 
general tendency for uKt to increase with pressure evident in our results is exactly the opposite of what was found 
in Ref. H^. In particular, they found a marked increase of aKr as p — > 0, which does not occur in our results. 

The thermodynamic Griineisen parameter 7 = V{dp/dE)v = uKtV / Cy plays an important role in high-pressure 
physics, because it relates the thermal pressure (i.e. the difference pth between p at given V and T and p at the 
same V but T — Q) and the thermal energy (difference i?th between E at given V and T and E at the same V but 
T = 0) . Assumptions about the value of 7 are frequently used in reducing shock data from Hugoniot to isotherm. If 
one assumes that 7 depends only on V , then the thermal pressure and energy are related by: 

pt^y = lEt^ , (23) 

a relation known as the Mie- Griineisen equation of state. At low temperatures, where only harmonic phonons 
contribute to i?th and pth, 7 should indeed be temperature independent above the Debye temperature, because 
i?th — ik^T per atom, and pthl^ = — 3fcBr(ilna'/(iln — ?ik^T^t\i, so that 7 = 7ph, which depends only on V . But 
in high-temperature Fe, the temperature independence of 7 will clearly fail, because of electronic excitations (and 
anharmonicity) . 

Our results for 7 (Fig. indicate that it varies rather little with either pressure or temperature in the region of 
interest. At temperatures below ca. 4000 K, it d ecrea ses with increasing pressure, as expected from the behaviour 



of the phonon Griineisen parameter 7ph (see Sec. IV B ). This is also expected from the often- used empirical rule of 
thumb 1^ 7 ~ {V/VqY, where Vq is a reference volume and g is a constant exponent usually taken to be roughly 
unity. Since V decreases by a factor of about 0.82 as p goes from 100 to 300 GPa, this empirical relation would make 
7 decrease by the same factor over this range, which is roughly what we see. However, the pressure dependence of 7 
is very much weakened as T increases, until at 6000 K 7 is almost constant. 

Our results agree quite well with those of Wasserman et al. in giving a value 7 ~ 1.5 at high pressures, although 
once again their calculations are limited to the low-temperature region T < 3000 K. But at low pressures there is a 
serious disagreement, since they find a strong increase of 7 to values of well over 2.0 as p — > 0, whereas our values 
never exceed 1.6. 



VII. DISCUSSION AND CONCLUSIONS 



Our primary interest in this work is in the properties of h.c.p. iron at high pressures and temperatures, but in 
order to investigate them using ab initio methods we have needed to make technical developments, which have a wider 
significance. The major technical achievement is that we have been able to calculate the ab initio free energy and other 
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thermodynamic properties with completely controlled statistical-mechanical errors, i.e. errors that can be reduced to 
any required extent. Anharmonicity and thermal electronic excitations are fully included. The attainment of high 
precision for the electronic and harmonic parts of the free energy has required no particular technical innovations, 
though careful attention to sources of error is essential. The main innovation is in the development of well optimised 
reference systems for use with thermodynamic integration in the calculation of the anharmonic part, without which 
adequate precision would be impossible. With the methods wc have developed, it becomes unnecessary to approximate 
the electronic structure with semi-empirical representations, or to resort to the statistical-mechanical approximations 
that have been used in the past. 

We have assessed in detail the precision achieved in the various parts of the free energy. There are two kinds of 
errors: those incurred in the calculation of the free energies themselves, and those produced by fitting the results to 
polynomials. We have seen that the errors in calculating the perfect-lattice free energy fperf are completely negligible, 
though there may be small fitting errors of perhaps 1 meV/atom. In the harmonic part Fharm, the calculational errors 
are ca. 3 meV/atom, most of which comes from spatial truncation of the force-constant matrix; the fitting error 
for -Fharm are of about the same size. The most serious errors are in the anharmonic part i^anhaim, and these are 
ca. 5 meV/atom in the calculation and ca. 4 meV/atom in the fitting. The overall technical errors therefore amount 
to ca. 15 meV/atom, which is slightly larger than our target of 10 meV/atom. 

We stress that the precision just quoted does not take into account errors incurred in the particular implementation 
of DFT (PAW in the present work), for example the error associated with the chosen split between valence and core 
states. Such errors can in principle be systematically reduced, but we have not attempted this here. Nor does it 
account for the inaccuracy of the chosen Ey^c, or for the neglect of the temperature dependence of E^c- We shall 
attempt to assess errors of this type in our separate paper on the melting properties of Fe. 

The most direct way to test the reliability of our methods is comparison with shock data for p{V) on the Hugo- 
niot [ pl| , so it is gratifying to find close agreement over the pressure range of interest. The closeness of this agreement 
is inherently limited by the known inaccuracies of the GGA employed, and we have shown that the discrepancies 
are of the expected size. An important prediction of the calculations is the temperature T{p) on the Hugoniot, since 
temperature is notoriously difficult to obtain in shock experiments. Our results support the reliability of the shock 
temperatures estimated by Brown and McQueen |^,^, and, in agreement with Wasserman et al. we find that 
the temperatures of Yoo et al. are too high by as much as 1000 K. This incidentally lends support to the reliability 
of the Brown and McQueen estimate of ca. 5500 K for the melting temperature of Fe at 243 GPa. The situation is 
not so satisfactory for the adiabatic bulk modulus Kg on the Hugoniot, since our ab initio values seem to be ca. 8 % 
above the shock values. But it should be remembered that even at ambient conditions ab initio and experimental 
bulk moduli frequently differ by this amount. The difficulties may be partly on the experimental side, since even for 
b.c.c. Fe at ambient conditions, experimental Ks values span a range of 8 %. 

Our calculations fully confirm the strong influence of electronic thermal excitations | |53| , |35[ |. At the temperatures 
T ^ 6000 K of interest for the Earth's core, their contribution to the specific heat is almost as large as that due to 
lattice vibrations, in line with previous estimates. They also have a significant effect on the Griineisen parameter 7, 
which plays a key role in the thermodynamics of the core, and is poorly constrained by experiment. Our finding that 
7 decreases with increasing pressure for T < 4000 K accords with an often- used rule of thumb ||5^ , but electronic 
excitations completely change this behaviour at core temperatures T ~ 6000 K, where 7 has almost constant values 
of ca. 1.45, in accord with experimental estimates in the range 1.1 to 1.6 [|l],|0). Compar ison with the earlier tight- 
binding calculations of Wasserman et al. both for 7 and for the quantity aKx is rather disquieting. Although 
a full comparison is hindered by the fact that they report results only for the low-temperature region T < 3000 K, 
we find two kinds of disagreement at low pressure. First, they find an increase of oKt as p ^ 0, whreas at low 
temperatures we find the opposite. Even more seriously, their strong increase of 7 as p ^ is completely absent 
in our results. Our calculations are more rigorous than theirs, since we completely avoid their statistical-mechanical 
approximations, as well as being fully self-consistent on the electronic-structure side. The suggestion must be that 
their approximations lead to significantly erroneous behaviour at low pressures. In pursuing this further, it would be 
very helpful to know what their methods predict in the high-temperature region relevant to the Earth's core. 

The present work forms part of a larger project on both pure Fe and its alloys with S, O, Si and H in the solid 
and liquid states. In a separate paper, we shall demonstrate that the thermodynamic integration technique employed 
here can also be used to obtain the fully ab initio free energy and other thermodynamic functions of liquid Fe over a 
wide range of states, with a precision equal to what has been achieved here for the solid. From the free energies of 
solid and liquid, we are then able to determine the ab initio melting curve and the entropy and volume of fusion as 
functions of pressure. 

In summary, we have presented extensive ab initio calculations of the free energy and a range of other thermodynamic 
properties of iron at high pressures and temperatures, in which all statistical-mechanical errors are fully under control. 
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and a high (and quantified) precision has been achieved. We find close agreement with the most rehable shock data. 
Ab initio values are provided for important, but experimentally poorly determined quantities, such as the Griineisen 
parameter. The free energy results provide part of the basis for the ab initio determination of the high-pressure 
melting properties of iron, to be reported elsewhere. 
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Vjk^) T(K) Fref - fharm (eV) JVib " J^ref (eV) Fyib " J'harm (eV) 

9.17 1000 4.448 -4.455 -0.007 

1500 4.437 -4.449 -0.012 

2500 4.416 -4.440 -0.024 

3000 4.405 -4.441 -0.036 

3500 4.403 -4.444 -0.041 

4000 4.402 -4.453 -0.051 

4500 4.385 -4.456 -0.071 

8.67 1000 4.961 -4.968 -0.007 

2000 4.933 -4.950 -0.017 

3000 4.913 -4.941 -0.028 

3500 4.906 -4.941 -0.035 

4000 4.902 -4.939 -0.042 

4500 4.893 -4.951 -0.058 

5000 4.888 -4.954 -0.067 

8.06 2250 5.683 -5.704 -0.021 

3400 5.655 -5.687 -0.032 

4500 5.634 -5.691 -0.057 

5000 5.626 -5.687 -0.061 

7.5 2500 6.530 -6.554 -0.024 

4000 6.490 -6.531 -0.041 

5000 6.468 -6.525 -0.057 

6000 6.448 -6.533 -0.085 

6.97 3000 7.531 -7.549 -0.018 

4500 7.483 -7.522 -0.039 

5500 7.460 -7.523 -0.063 

6000 7.449 -7.514 -0.065 

6500 7.438 -7.519 -0.081 

7000 7.428 -7.523 -0.095 



TABLE I. Anharmonic contribution Fanharm = -FVib — Fharm to the ab initio free energy of h.c.p. Fe at a set of atomic 
volumes V and temperatures T, calculated as a sum of the free energy difference F^gf — Fharm between reference and ab initio 
systems and the difference Fvib — F-ef between full ab initio and reference systems. All quantities are per atom. 
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FIG. 1. Electronic specific heat per atom Cperf of the rigid perfect lattice of h.c.p. Fe (units of Boltzmann's constant /cb) as 
a function of temperature for atomic volumes: 7.0 ( ), 8.0 ( ), 9.0 A^ (- - -) and 10.0 A^ (■■■). 

FIG. 2. Electronic density of states of h.c.p. Fe calculated at atomic volumes of 7.0 A'^ ( ) and 10.0 A^(- ■ ■). Energy is 

referred to the Fermi energy Ef . 

FIG. 3. Electronic thermal pressure Appcrf of the rigid perfect lattice of h.c.p. Fe as a function of atomic volume V for 
T = 2000 ( ), 4000 ( ) and 6000 K (- - -). 

FIG. 4. Phonon dispersion relations of h.c.p. Fe calculated at atomic volumes V = 8.67 (left panel) and 6.97 A'^ (right 
panel). Frequencies calculated directly from DFT at the two volumes are shown as solid curves. In left panel, dashed curves 
give frequencies from empirical inverse-power model (see text). In right panel, dashed curves show DFT frequencies for 
V — 8.67 A^ graphed in left panel but scaled by the factor 1.409. 

FIG. 5. Geometric-mean phonon frequency dj of h.c.p. Fe as a function of atomic volume V for T — 2000 ( ), 4000 ( ) 

and 6000 K ( — ). The natural logarithm of the two quantities is plotted, with uj in units of rad s~^ and V in units of A"^. 

FIG. 6. The harmonic thermal pressure Pharm as a function of atomic volume V for T = 2000 ( ), 4000 ( ) and 6000 K 

(---). 

FIG. 7. Experimental and ab initio Hugoniot pressure p as a function of atomic volume V. Symbols show the measurements 
of Brown and McQueen |Q. Solid curve is ab initio pressure obtained when calculated equilibrium volume of b.c.c. Fe is used 
in the Hugoniot-Rankine equation; dotted curve is the same, but with experimental equilibrium volume of b.c.c. Fe. 

FIG. 8. Experimental and ab initio temperature as a function of pressure on the Hugoniot. Black circles with error bars and 
white diamonds are measurements due to Yoo et al. and Brown and McQueen respectively. Solid and dashed curves 
are ab initio results obtained using theoretical and experimental b.c.c. volumes. 

FIG. 9. Experimental and ab initio adiabatic bulk modulus Ks on the Hugoniot. Diamonds and pluses are measurements 
due to Jeanloz |^ and Brown and McQueen jil] respectively. Solid and dashed curves are ab initio results obtained using 
theoretical and experimental b.c.c. volumes. 

FIG. 10. Experimental and ab initio thermal expansivity on the Hugoniot. Diamonds are measurements due to Jeanloz [ po[ . 
Solid and dashed curves are ab initio results obtained using theoretical and experimental b.c.c. volumes. 

FIG. 11. Total constant-volume specific heat per atom Cv (units of fee) of h.c.p Fe as a function of pressure from present ab 
initio calculations at T = 2000 ( ), 4000 ( ) and 6000 K (- - -). 

FIG. 12. Ab initio thermal expansivity a as a function of pressure on isotherms T = 2000 ( ), 4000 ( ) and 6000 K (- - 

-). Black circle with error bar is experimental value of Duffy and Ahrens at T = 5200 ± 500 K. Diamonds are DAC values 
due to Boehler fin] for temperatures between 1500 and 2000 K. 



FIG. 13. Ab initio values of product of thermal expansion coefficient ct and isothermal bulk modulus Kt as a function of 
pressure P at T = 2000 ( ), 4000 ( ) and 6000 K (- - -). 



FIG. 14. Ab initio Griineisen parameter 7 as a function of pressure at T = 2000 ( ), 4000 ( ) and 6000 K ( — ). 



18 



FIGURE 1 




T(K) 



19 



FIGURE 2 

1 1 1 1 1 r 




E - EF(eV) 



20 



Thermal Pressure (GPa) 



O 



cn 




22 




23 



FIGURE 6 




24 




25 



FIGURE 8 
9000 I \ \ 




P (GPa) 



26 




27 




28 




29 



FIGURE 12 




50 100 150 200 250 300 350 



P (GPa) 



30 



FIGURE 13 

1 1 1 r 




I I I I I I I 

50 100 150 200 250 300 350 

P (GPa) 



31 



FIGURE 14 




50 100 150 200 250 300 350 

P (GPa) 



32 



